AD-A114  Oil 

UNCLASSIFIED 


ROYAL  AIRCRAFT  ESTABLISHMENT  FARNBOROUGH  (ENGLAND)  F/G  S/2 

THE  USE  OF  TERRAIN  HEIGHT  INFORMATION  FOR  IMPROVING  THE  ACCURAC—ETC ( U> 
JAN  82  P  A  ROBERTS 
RAE-TM-SPACE-297 


1 


DRIC-BR-82645 


NL 


I 


ROYAL  AIRCRAFT  ESTABLISHMENT 


THE  USE  OF  TERRAIN  HEIGHT  INFORMATION  FOR  IMPROVING  THE 
ACCURACY  OF  CLASSIFICATION  OF  LANDSAT  DATA 

by 

P .  A .  Roberts 


January  1982 


This  doeuamt  has  Gm  appwwd 
lor  public  roloaao  and  scdoi  ttt 
diotrlbudoa  it  unUudtcd. 


DTIC 

EUECTE 
APR  2  9 1982 

E 


04  27  .  027 


82 


ROYAL  AIRCRAFT  ESTABLISHMENT 
Technical  Memorandum  Space  297 
Received  for  printing  12  January  1982 

THE  USE  OF  TERRAIN  HEIGHT  INFORMATION  FOR  IMPROVING  THE 
ACCURACY  OF  CLASSIFICATION  OF  LANDS AT  DATA 

by 

P.  A.  Roberts 

\  SUMMARY 

V  One  potentially  important  area  of  work  in  remote  sensing  is  the  use 
of  topographic  data  in  conjunction  with  satellite  imagery  to  improve  the 
accuracy  of  classification  algorithms.  This  Memorandum  is  concerned  with 
techniques  for  using  terrain  height  data  in  conjunction  with  Landsat  imagery. 
An  algorithm  is  described  that  transforms  digitised  height  contours  into  an 
equispaced  grid  of  points.  A  method  is  also  described  which  utilises 
this  height  matrix  to  derive  the  magnitude  and  direction  of  the  slope  of 
the  terrain,  information  which  can  in  turn  be  used  to  calculate  the  varia¬ 
tion  in  direct  solar  illumination  over  an  area.  Finally,  techniques  are 
described  for  using  terrain  data  to  improve  the  accuracy  of  multispectral 
classification  of  Landsat  data. ^ 
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I  INTRODUCTION 

Much  more  information  can  often  be  obtained  from  remotely  sensed  data  if  it  is 
used  in  association  with  existing  ground  truth  information.  This  paper  is  concerned 
with  techniques  for  using  topographic  information  from  maps,  in  conjunction  with  imagery 
of  land  surfaces.  An  algorithm  is  described  that  can  be  used  to  transform  digitised 
contour  height  information  into  an  equispaced  grid  of  points,  a  format  that  is  suitable 
for  easy  manipulation  and  processing  and  for  comparison  with  many  types  of  digital 
imagery.  The  height  matrix  can  in  turn  be  used  to  derive  the  magnitude  and  direction 
of  the  slope  of  the  terrain  and  the  variation  in  the  direct  solar  illumination  arising 
from  it.  These  quantities  can  all  be  displayed  as  grey  tone  images  for  independent 
interpretation  or  interpretation  in  conjunction  with  other  imagery.  In  addition,  the 
intensity  of  Landsat  pixels  can  be  adjusted  so  as  to  compensate  for  variations  in  direct 
solar  illumination  arising  from  undulating  terrain. 

In  this  Memorandum  the  factors  that  need  to  be  taken  into  account  when  digitising 
a  contour  map  are  first  briefly  discussed.  This  is  then  followed  by  a  description  of 
a  technique  for  producing  a  terrain  height  matrix  from  a  contour  height  map.  Subsequent 
sections  deal  with  the  parameters  that  can  be  obtained  from  the  terrain  height  matrix 
and  the  display  of  these  quantities  as  grey  tone  images.  The  final  part  of  this 
Memorandum  describes  how  terrain  height  data  can  be  used  to  try  to  improve  the  accuracy 
of  multispectral  classification  of  Landsat  data.  In  order  to  illustrate  these 
techniques  a  small  test  area  in  Southern  England  has  been  selected. 


2  DIGITISATION  OF  CONTOURS 

Often  the  investigator  will  not  have  the  contour  information  already  digitised, 
but  rather  will  only  have  a  paper  copy  of  the  map  thus  necessitating  digitisation  of  the 
contours.  It  is  necessary  therefore  for  the  investigator  to  decide  at  what  interval 
around  contours  points  are  to  be  digitised.  The  evaluation  of  this  interval  involves 
an  estimate  being  made  of  the  interval  along  a  contour  that  can  be  represented  by  a 
straight  line.  This  interval  is  dependent  upon  the  grid  spacing  of  the  terrain  height 
matrix  and  the  curvature  of  the  contour  in  question;  the  smaller  the  grid  spacing  or  the 
greater  the  curvature  of  the  contour,  the  smaller  this  interval  must  be  in  order  to 
represent  the  contour  to  the  required  accuracy.  More  precisely,  if  d  is  the  digitisa¬ 
tion  interval,  r  is  the  radius  of  curvature  of  the  contour  (which  generally  is  a  func¬ 
tion  of  position  along  the  contour)  and  s  is  the  maximum  permissible  deviation  of  a 
straight  line,  joining  two  points  on  the  contour,  from  the  contour,  then  d  ^  /8rs  , 
where  d  is  such  that  the  angle  the  contour  turns  through  is  small.  The  author  has 
found  that  for  his  work  it  is  adequate  that  s  be  half  the  grid  spacing,  but  this  may 
not  be  acceptable  for  all  applications. 

If  the  digitisation  interval  is  greater  than  d  then  the  contour  is  under¬ 
digitised,  and  if  it  is  less  it  is  over-digitised.  The  only  consequence  of  the  latter 
is  that  the  number  of  data  points  describing  the  contour  is  greater  than  necessary. 
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Contour  fitting  algorithms  that  are  considerably  better  than  linear  interpolation  can  be 
employed  and  should  be  adopted  when  contours  are  under-digitised.  A  description  of  one 
such  algorithm  can  be  found  in  Ref  I . 

The  digitisation  of  contours  can  take  either  of  the  following  forms  (amongst 
others) : 

(i)  Constant  separation  between  data  points. 

(ii)  Variable  density  system  with  the  interval  between  data  points  depending  upon 
the  curvature  of  the  contour.  This  reduces  the  number  of  points  required  to 
describe  a  contour. 

If  it  is  decided  to  opt  for  a  constant  separation  between  data  points  then  the  map 
will  be  correctly  digitised  if  the  evaluation  of  the  digitisation  interval  is  based  upon 
the  smallest  radius  of  curvature.  The  variable  density  method  requires  a  continual 
evaluation  of  the  curvature  of  the  contour  being  digitised  and  is  best  suited  to  auto¬ 
mated  digitisation.  If  the  digitisation  is  being  performed  manually  using  this  method 
then  the  investigator  must  make  a  subjective  estimate  of  the  required  digitisation 
interval . 

3  TRANSFORMING  DIGITISED  CONTOURS  INTO  A  HEIGHT  MATRIX 

The  basic  problem  in  generating  a  height  matrix  is  to  fit  a  surface  to  the  contour 
data  such  that  every  contour  lies  on  the  surface.  Since  it  is  assumed  that  there  will 
be  no  further  information  available  as  to  the  nature  of  the  surface,  the  intuitively 
desirable  conditions  of  smoothness  and  continuity  are  adopted.  This  surface  can  then  be 
sampled  at  the  appropriate  points  to  produce  the  required  grid.  It  should  be  remembered 
that  the  surface  has  to  be  constrained  such  that  between  contours  it  remains  within  the 
range  of  values  defined  by  the  contours. 

The  method  adopted  does  not  perform  a  true  two-dimensional  fit  but  instead  takes 
the  average  of  two  orthogonal  one  dimensional  fits.  This  involves  two  orthogonal  sets 
of  cuts  being  taken  across  the  contour  map;  the  points  of  intersection  between  these 
grid  lines  define  the  points  at  which  interpolated  data  is  required.  For  each  grid  line 
the  points  of  intersection  of  the  line  with  the  contours  are  calculated  and  then  inter¬ 
polation  is  carried  out  along  the  grid  line  to  obtain  values  at  the  grid  points  lying  on 
the  line.  Also  produced  are  weighting  values  reflecting  the  confidence  in  the  inter¬ 
polated  values.  This  procedure  results  in  two  grids  of  points  being  obtained  correspond¬ 
ing  to  the  two  sets  of  cuts,  together  with  two  grids  containing  the  associated  weighting 
values.  These  weighting  values  are  used  to  combine  the  grids  to  produce  the  final  height 
matrix. 

It  is  necessary  to  take  two  mutually  orthogonal  sets  of  cuts  across  the  map  since 
it  may  be  that  some  contours  run  parallel  to  one  of  the  sets  of  cuts  and  thus  would  be 
ignored  by  that  set  unless  they  fell  exactly  on  the  cuts.  Contours  missed  by  one  set 
of  cuts  are  therefore  incorporated  in  the  orthogonal  set.  In  this  way  the  orientation 
of  the  cuts  with  respect  to  the  contour  map  is  not  important,  whilst  it  may  be  crucial 
if  only  one  set  of  cuts  were  to  be  taken. 
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When  converting  a  contour  map  into  a  height  matrix  it  is  necessary  to  define  a 
grid  spacing.  It  should  be  remembered  that  the  smaller  the  grid  spacing  the  more 
accurately  the  contours  can  be  reconstructed  from  the  height  matrix,  and  in  any  case  it 
should  be  no  greater  than  half  the  scale  size  of  the  smallest  feature  it  is  required  to 
preserve  in  the  transformation.  Fig  1  shows  a  5  km  square  area  of  Hampshire  after 
digitisation  of  the  contours  and  spot  heights  from  a  1:25000  Ordnance  Survey  Map  (sheet 
SU54).  This  map  has  been  converted  into  a  grid  500  *  500  using  a  10  m  grid  spacing. 

The  accuracy  of  the  height  matrix  can  be  ascertained  by  comparison  with  the  spot  heights 
for  the  area,  although  it  should  be  realised  that  beyond  the  outermost  contours  the 
interpolation  technique  is  not  very  reliable.  Values  can  be  derived  from  the  height 
matrix  at  the  positions  of  the  spot  heights  by  using  the  16  surrounding  matrix  values. 

A  comparison  between  the  two  sets  of  values  is  given  in  Table  1 .  As  is  to  be  expected, 
there  are  no  differences  greater  than  the  contour  interval  of  25  ft.  The  mean  differ¬ 
ence  between  the  two  sets  of  heights  is  only  0.04,  a  value  that  is  probably  rather 
fortuitous  for  a  sample  of  only  21  points.  The  standard  deviation  of  the  differences 
is  7.1  ft. 


A  much  more  detailed  description  of  the  algorithm  used  to  generate  the  height 
matrix  from  contours  is  given  in  Ref  1.  This  reference  also  includes  a  description  and 
a  listing  of  a  computer  program  that  can  be  used  for  evaluation  of  the  algorithm. 
Further  information  and  refinements  of  the  algorithm  are  to  be  found  in  Ref  2. 


4  QUANTITIES  DERIVED  FROM  THE  HEIGHT  MATRIX 

Using  the  height  matrix  it  is  possible  to  derive  the  magnitude  and  direction  of 
the  gradient  of  the  terrain.  Henceforth  in  this  Memorandum  the  magnitude  of  the  gradient 
will  be  referred  to  as  slope  and  the  direction  as  aspect.  The  slope  and  aspect  are 
derived  from  the  x  and  y  direction  gradients  ra^  and  m^  that  are  obtained  from  a 
least  squares  fit  of  a  two  dimensional  second  order  polynomial  over  a  small  n  x  n 
matrix  of  points.  This  least  squares  fitting  is  considerably  simplified  by  the  fact 
that  the  points  are  an  equispaced  grid  and  the  point  of  interest  is  the  one  at  the 
centre.  The  least  squares  fit  is  performed  using  either  3x3#5x5or7x7  points 
depending  upon  how  much  smoothing  is  considered  desirable. 


The  slope  and  aspect  are  obtained  from  the  x  and  y  direction  gradients  m 


and  m  as  follows: 

y 


-/  n  r\ 

slope  *  tan  Uta^  x  I 


-i/V 

aspect  *  tan  I  —  ] 

my  / 


undefined 


if  slope  ^  0 


if  slope  »  o. 


The  direct  illumination  from  the  Sun  at  a  point  on  the  ground  can  in  turn  be 
obtained  from  the  slope  and  aspect,  viz: 


L 
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Ig  «  {cos (a  -  y)  cos  B  sin  <5  +  sin  B  cos  5)1^ 

where  a  is  the  Sun  azimuth,  B  is  the  Sun  elevation,  y  the  aspect,  6  the  slope 
and  Iq  the  solar  illumination  on  a  unit  area  perpendicular  to  the  direction  of  the 
illumination;  a  and  y  are  measured  from  north.  It  must  be  emphasised  that  this 
equation  only  refers  to  direct  illumination  and  neglects  totally  the  not  insignificant 
contribution  that  will  arise  from  backscattered  radiation.  A  quantity  that  is  particu¬ 
larly  useful  in  investigating  the  variation  in  brightness  arising  from  the  surface 
topography  is  given  by: 


n  sin  6  cos (a  -  y)  .  x 

d  tan  B 

This  quantity,  which  henceforth  will  be  referred  to  as  the  variation  in  direct 
solar  illumination,  is  unity  for  flat  terrain,  greater  than  unity  for  terrain  inclined 
towards  the  Sun,  less  than  unity  for  terrain  inclined  away  from  the  Sun  and  negative  for 
terrain  in  shadow.  Also,  the  variation  in  over  a  given  area  increases  as  the  Sun 

elevation  decreases. 

Neither  of  these  equations  takes  into  account  the  possibility  that  some  parts  of 
the  terrain  may  be  shaded  from  direct  illumination  by  others.  Thus,  these  equations  are 
only  valid  for  all  values  of  the  terrain  aspect  and  Sun  azimuth  angles  when  the  Sun 
elevation  angle  is  greater  than  the  steepest  slope,  ie  B  >  <5  . 

5  DISPLAYING  TERRAIN  PARAMETERS  AS  GREY  TONE  IMAGES 

Height,  slope  and  variation  in  direct  solar  illumination  present  no  problems  in 
being  displayed  as  grey  tone  images,  since  lower  values  can  be  assigned  to  lower  pixel 
intensities  and  higher  values  to  higher  pixel  intensities.  Aspect  can  also  be  displayed 
as  a  grey  tone  image  although  the  interpretation  of  such  an  image  is  rather  more 
difficult.  This  is  because: 

(1)  The  aspect  can  vary  between  0°  and  360°,  however  0°  and  360°  represent  the 
same  direction. 

(?)  The  aspect  for  flat  terrain  is  undefined. 

(3)  The  value  of  the  aspect  becomes  more  unreliable  as  the  slope  becomes  smaller. 

One  convention  that  can  be  adopted  for  displaying  aspect  and  the  one  that  is 
adopted  here,  although  there  is  some  ambiguity  in  the  interpretation  of  the  minimum  and 
maximum  permissible  image  intensities,  is  to  assign  an  aspect  value  of  0°  to  an  intensity 
of  zero,  360°  to  the  maximum  image  intensity,  and  flat  terrain  to  a  zero  image  intensity. 
Other  conventions  can  of  course  be  envisaged. 

When  considering  the  height  matrix  as  a  grey  tone  image  it  must  be  remembered  that 
each  value  is  associated  with  the  centre  of  a  pixel  and  does  not  represent  the  average 
value  over  a  pixel.  Where  the  terrain  is  changing  slowly  compared  with  the  size  of  a 
pixel  then  the  values  of  the  height  matrix  will  be  a  good  approximation  to  the  average 
value  over  a  pixel.  However,  where  the  terrain  is  changing  rapidly  compared  with  the 
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pixel  size  the  height  estimate  will  not  necessarily  be  a  good  approximation  to  the  mean 
value  over  a  pixel.  In  this  situation  the  height  matrix  must  first  be  evaluated  using 
a  grid  spacing  for  which  the  terrain  is  changing  slowly  and  then  a  suitable  average  taken 
to  obtain  an  estimate  of  the  mean  value  for  the  required  pixel  size.  This  technique  can 
also  be  applied  to  the  variation  in  direct  solar  illumination  but  it  cannot  be  applied  to 
the  slope  and  aspect.  For  these  quantities  it  is  necessary  to  evaluate  the  mean  x  and 
y  direction  gradients  in  the  manner  described  and  thence  to  derive  the  mean  slope  and 
aspect . 

As  an  example,  the  central  4  km  square  area  of  the  map  shown  in  Fig  1  has  been 
converted  to  a  height  matrix  using  a  grid  spacing  of  10  m.  This  was  then  used  to  derive 
the  height,  slope,  aspect  and  variation  in  direct  solar  illumination  corresponding  to 
50  m  pixels.  (The  reason  for  choosing  this  pixel  size  will  become  evident  in  the  next 
section.)  The  values  associated  with  the  pixels  have  been  obtained  by  using  the  25 
values  of  the  matrix  covering  each  pixel.  The  four  terrain  parameters  are  displayed  as 
grey  tone  images  in  Fig  2. 

6  APPLICATION  TO  LANDSAT  IMAGERY 

6 . 1  Methods  of  employing  terrain  data 

Recently,  an  increasing  amount  of  interest  has  been  given  to  the  possibility  of 
improving  the  accuracy  of  multispectral  classification  with  Landsat  imagery  by  taking 
into  account  effects  arising  from  surface  topography.  These  effects  arise  from  the  fact 
that  although  two  areas  may  have  the  same  cover  type,  and  thus  for  a  given  spectral  band 
would  be  expected  to  have  the  same  intensity,  in  practice  they  will  have  different 
intensities  if  the  terrain  is  such  that  the  two  areas  receive  differing  amounts  of 
illumination.  If  data  describing  the  surface  topography  is  available  it  should  be 
possible  to  reduce  this  topographic  effect  and  thus  improve  the  accuracy  of  classifica¬ 
tions.  This  can  be  achieved  by  two  methods:  one  method  involves  using  the  terrain  para¬ 
meters  as  extra  image  planes  and  the  other  involves  modifying  the  pixel  intensities  in 
the  Landsat  image. 

In  the  first  method,  the  four  terrain  parameters  can  be  considered  as  four 

additional  image  planes  to  be  used  in  conjunction  with  the  four  image  planes  associated 

with  the  four  spectral  bands  of  a  Landsat  scene.  All  eight  image  planes  can  be  analysed 

3 

together  either  on  an  interactive  image  processing  system  such  as  the  IDP3000  or  on  a 
computer.  Thus,  in  classification,  it  would  be  possible  to  consider  only  areas  exceed¬ 
ing  a  specified  height  or  only  flat  areas.  When  classifying  areas  of  similar  cover  type 
the  reliability  of  the  classification  may  be  improved  by  performing  a  series  of  classi¬ 
fications  each  involving  pixels  which  represent  areas  which  received  similar  amounts  of 
illumination.  This  can  be  achieved  by  using  the  image  plane  containing  the  variation  in 
direct  solar  illumination  in  conjunction  with  the  four  Landsat  image  planes. 

The  second  method  of  improving  the  reliability  of  classif ication  is  rather  more 
difficult  and  involves  more  assunptions.  Here,  corrections  are  made  to  the  intensity 
of  pixels  so  as  to  account  for  both  the  variation  in  illumination  and  for  the  scattering 
properties  of  the  surface.  The  total  illumination  consists  of  direct  solar  illumination 
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and  a  contribution  from  a  diffuse  component  comprised  primarily  of  scattered  skylight, 
although  there  is  a  small  contribution  from  light  scattered  onto  the  surface  by  surround- 
ing  ground*  The  diffuse  component  has  been  shown  to  be  a  small  proportion  of  the  total 
illumination,  ie  less  than  12%  of  the  tctal  for  each  of  the  four  Landsat  bands  for  a  Sun 
elevation  of  34°  under  clear  sky  non-hazy  conditions  (Ref  4),  conditions  which  will  be 
satisfied  by  many  Landsat  scenes  that  would  generally  be  considered  suitable  for  classify¬ 
ing  surface  features.  A  further  complication  arises  from  the  fact  that  most  natural 
surfaces  have  preferred  directions  of  scattering  and  this  should  be  taken  into  account. 
This  is  difficult  because  the  scattering  properties  of  surfaces  are  dependent  upon  their 
composition  and  roughness  and  hence  are  hard  to  model.  Further,  the  nature  of  the 
scattering  may  be  heavily  dependent  on  transient  conditions,  eg  the  presence  of  surface 
water  on  vegetation  after  rain  or  a  heavy  dew.  Because  of  the  complexity  of  the  problem 
the  far  from  insigificant  effect  due  to  the  scattering  properties  of  the  surface  are 
ignored  here.  A  detailed  discussion  of  the  nature  of  scattering  by  natural  surfaces 
can  be  found  in  Ref  5.  Thus,  assuming  that  all  surfaces  scatter  light  equally  in  all 
directions,  pixel  values  can  be  adjusted  so  as  to  take  account  of  variations  in  the 
surface  illumination  by  using  the  formula: 


where  1^  =  uncorrected  pixel  intensity 

Ic  =  corrected  pixel  intensity 

and  r  =  diffuse  component  expressed  as  a  fraction  of  the  direct  illumination  on  a 

horizontal  surface 

It  must  be  realised  that  attempting  to  correct  pixels  in  shadow  will  not  be  expected  to 
produce  much  improvement  in  classification  because  of  the  small  range  of  pixel  values  in 
such  areas. 

To  illustrate  the  effects  of  this  technique,  a  Landsat  scene  covering  the  area 
illustrated  in  Fig  I  has  been  selected.  Each  of  the  four  Landsat  bands  has  been 
destriped,  geometrically  corrected  and  brought  into  registration  with  the  National  Grid. 

A  subscene  corresponding  to  the  central  4  km  square  of  Fig  1  has  been  extracted  from  the 
main  scene  to  enable  direct  comparison  to  be  made  with  the  terrain  parameters  for  the 
area.  The  pixel  size  for  the  Landsat  data  has  been  selected  to  be  50  ra  square  so  as  to 
preserve  as  much  as  possible  the  accuracy  of  the  original  data.  (The  Landsat  image  has 
to  be  resampled  during  the  geometrical  correction  process.  For  further  details  see  for 
example  Ref  6.)  A  late  Autumn  scene  (25  November  1 973)  with  a  Sun  elevation  of  only  15°  i 

has  been  selected  so  as  to  highlight  the  effect  of  the  variation  in  direct  solar  illumi¬ 
nation  over  undulating  terrain.  In  comparison,  all  but  10%  of  the  pixels  are  associated 
with  slopes  less  than  5°  and  the  mean  slope  is  only  2.7°.  Each  Landsat  pixel  was  modi¬ 
fied  using  the  above  equation  with  a  value  of  0.14  being  adopted  for  r  .  Fig  3  shows 
the  Landsat  band  7  image,  corresponding  to  the  central  4  kra  square  of  Fig  1,  both  before 
and  after  correction.  Fig  3a  is  the  original  data  after  a  suitable  contrast  stretch 
has  been  performed  and  Fig  3b  is  the  same  image  but  after  the  topographic  correction 
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has  been  applied.  As  can  be  seen  from  a  comparison  of  the  two  images,  the  correction 
has  quite  a  marked  effect,  which  is  perhaps  especially  important  since  the  test  area  is 
typical  of  many  areas  of  gently  undulating  terrain  found  in  Southern  England.  The  effect 
will  be  expected  to  be  even  greater  in  areas  where  steeper  gradients  exist  and  farther 
north  where  lower  Sun  elevations  are  encountered. 

In  the  above  analysis  the  Sun  elevation  angle  has  been  assumed  to  be  the  same  over 
all  the  test  area.  This  is  a  reasonable  assumption  here  since  the  test  area  is  only 
4  km  ^  4  km,  but  allowance  should  be  made  for  the  variation  in  Sun  elevation  angle  if 
large  areas  are  considered.  For  example,  over  a  Landsat  scene  of  the  UK  the  variation 
is  of  the  order  of  2°,  which  can  give  rise  to  significant  effects  for  low  Sun  angles. 

6.2  Comparison  of  different  classification  techniques 

The  extent  of  any  improvement  in  the  accuracy  of  classification  using  Landsat  data 
when  combined  with  terrain  data  can  only  be  ascertained  by  comparison  with  classifica- 
tions  which  do  not  involve  the  use  of  terrain  data,  and  by  comparison  of  these  classi¬ 
fications  with  ground  truth  data.  Before  describing  the  results  of  different  classifica¬ 
tions  it  is  important  to  understand  the  classification  procedure  that  is  employed.  First,  an 
area  is  selected  (of  ten  referred  to  as  the  1  training*  area)  over  which  it  is  known  that  the 
cover  type  of  interest  lies  and  it  is  assumed  that  all  pixels  within  this  area  are  part 
of  the  class  of  interest.  The  minimum  and  maximum  intensities  for  each  band  are  then 
determined  for  these  pixels  and  then  all  other  pixels  in  the  scene  are  examined  to 
determine  whether  or  not  they  lie  within  the  above  limits  and  if  they  do  they  are  then 
included  in  the  classified  area. 

It  was  realised  at  the  outset  that  the  test  area  was  too  small  and  the  terrain 
insufficiently  rugged  for  a  reliable  comparison  to  be  made  of  different  classification 
techniques.  The  test  area  was  initially  selected  because  it  would  not  involve  too  much 
effort  to  digitise  as  its  main  purpose  was  to  assist  in  developing  the  algorithms  for 
utilizing  terrain  data.  Problems  were  also  encountered  in  trying  to  establish  ground 
truth  to  verify  the  classifications.  It  was  decided  to  try  to  classify  areas  of  woodland 
since  reliable  information  as  to  the  extent  of  woodland  in  the  area  was  available  from 
the  Forestry  Commission.  Most  of  the  remainder  of  the  test  area  was  farmland.  An  out¬ 
line  was  therefore  drawn  around  the  main  areas  of  woodland  but  approximately  50  m  inside 
the  boundary  to  remove  uncertainties  that  might  arise  from  a  possible  misregistration 
between  the  Landsat  image  and  the  map  data.  (This  is  discussed  in  more  detail  later.) 
Using  the  outlines  drawn,  the  Landsat  pixels  that  were  associated  with  woodland  were 
determined.  A  further  difficulty  arose  because  the  woodland  was  largely  a  mixture  of 
coniferous  and  deciduous  trees  and  this  mixture  was  different  for  different  areas.  Also, 
as  the  Landsat  data  was  obtained  in  late  Autumn  the  different  species  of  deciduous  trees 
were  at  different  stages  of  shedding  their  foilage.  Because  of  all  these  difficulties 
it  was  realised  that  little  faith  could  be  put  in  any  conclusions  drawn  from  the  results 
obtained  but  nevertheless  the  exercise  in  itself  was  considered  to  be  useful. 

It  was  decided  to  compare  the  classifications  obtained  when  using  the  following 
three  types  of  data: 
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(1)  The  original  Landsat  data. 

(2)  The  original  Landsat  data  after  ratioing.  (In  this  case  each  band  is  divided 
by  the  sum  of  all  four  bands,  the  rationale  for  this  being  that  the  resultant 
four  image  planes  should  be  independent  of  the  level  of  illumination.  This 
method  was  chosen  for  examination  because  it  is  the  one  commonly  employed  for 
attempting  to  remove  the  topographic  effect.) 

(3)  The  illumination  corrected  Landsat  data. 

The  results  of  the  classification  are  given  in  Table  2.  The  second  column 
represents  the  total  number  of  pixels  in  the  subscene  that  are  classified  as  being  mixed 
woodland  whereas  the  third  column  represents  the  number  of  classified  pixels  coincident 
with  pixels  that  were  definitely  known  to  be  associated  with  mixed  woodland.  The  last 
column  represents  the  proportion  of  the  test  areas  of  mixed  woodland  correctly  classified. 
The  two  important  points  to  note  are  that  the  classifications  with  the  original  data  and 
the  illumination  corrected  data  are  equally  successful  and  that  the  ratio  method 
produced  much  poorer  results.  It  is  often  found  that  the  ratio  method  gives  inferior 
results  compared  with  classifications  using  the  original  data  because  although  in 
principle  it  removes  the  effect  of  different  levels  of  illumination  it  is  at  the  expense 
of  losing  intensity  information  which  might  be  essential  for  a  successful  classification. 

A  comparison  was  also  made  between  the  intensity  distributions  of  the  pixels  cover¬ 
ing  the  woodland  areas,  both  before  and  after  the  topographic  correction  was  applied. 

What  one  would  expect  to  find  if  the  correction  improves  the  data  for  classification 
is  that  the  variance  of  the  pixel  intensities  is  smaller  after  the  correction  is  applied. 
The  results  obtained  can  be  seen  in  Table  3.  These  show  that  after  the  correction  the 
variance  is  smaller  for  band  7,  unchanged  for  band  6,  greater  for  band  5  and  greater 
still  for  band  4.  The  reasons  for  this  are  not  clear.  One  possible  explanation  is  that 
the  assumption  concerning  the  relative  contribution  of  the  diffuse  skylight  is  incorrect. 
The  assumption  made  was  that  for  each  band  it  represented  12%  of  the  total  illumination 
on  flat  terrain.  However,  it  would  be  expected  that  the  relative  contribution  would  be 
greater  at  the  shorter  wavelengths  because  the  effect  of  scattering  by  dust  particles 
and  water  droplets  in  the  atmosphere  is  greater  at  these  wavelengths.  Also,  no  account 
has  been  taken  of  the  fact  that  the  diffuse  contribution  is  anisotropic,  this  being 
manifest  as  an  increase  in  intensity  near  the  horizon  and  in  the  circumsolar  region  of 
the  sky.  Another  possibility  is  that  account  needs  to  be  taken  of  the  scattering 
properties  of  the  type  of  woodland  under  consideration,  data  which  at  present  is  not 
available  to  the  author. 

6 . 3  Discussion 

Having  developed  the  relevant  techniques  for  using  terrain  data  with  Landsat 
imagery  the  next  step  is  to  carry  out  a  comprehensive  analysis  of  improvements  in  the 
accuracy  of  classification  arising  from  knowledge  of  the  terrain.  When  selecting  a  new 
test  area  for  this  analysis  four  criteria  will  need  to  be  satisfied: 


SI*  297 


1  1 


(1)  The  area  should  be  much  larger  than  the  test  area  used  here. 

(2)  The  terrain  should  be  considerably  more  rugged  than  for  the  test  area  used 
here . 

(3)  Very  accurate  ground  truth  information  should  be  available  to  verify 
classifications. 

(4)  Very  good  terrain  data  is  required,  eg  height  contours  every  25  ft.  It  would 
also  be  interesting  to  determine  how  any  improvement  in  classification  changed 
as  the  quality  of  the  map  data  used  to  generate  the  height  matrix  was  reduced. 

At  present,  the  most  promising  method  for  employing  terrain  data  to  improve 
classification  accuracy  is  for  the  variation  in  direct  illumination  to  be  used  as  an 
extra  image  plane.  Trying  to  correct  an  image  so  that  it  appears  as  it  would  were  the 
terrain  flat  rather  than  undulating  has  the  disadvantage  that  both  the  contribution  of 
the  diffuse  skylight  and  the  scattering  properties  of  the  cover  type  of  interest  need  to 
be  known,  two  quantities  that  are  rather  awkward  to  determine. 

One  problem  that  has  been  ignored  so  far  is  that  arising  from  possible  mis¬ 
registration  between  the  Landsat  data  and  the  map  data  resulting  in  a  height  (or  other 
terrain  parameter)  estimate  associated  with  a  particular  pixel  not  being  associated  with 
that  pixel  but  with  a  nearby  pixel  instead.  The  misregis tration  arises  from  the 
positional  uncertainty  in  both  the  Landsat  data  and  the  map  data.  The  accuracy  of  the 
Landsat  data  after  geometrical  correction  is  typically  about  50  m,  this  value  represent¬ 
ing  the  root  mean  square  value  of  the  differences  between  the  positions  of  ground  control 
points  on  the  map  and  on  the  Landsat  image.  For  an  Ordnance  Survey  map,  the  vertical 
accuracy  of  the  contours  can  be  taken  as  being  approximate ly  equivalent  to  a  quarter  contour 
interval,  however  Ref  7  should  be  consulted  for  a  precise  definition  of  contour  accuracy. 
Thus,  the  misregistration  is  likely  to  be  greatest  for  flat  and  undulating  terrain  where 
contours  are  widely  spaced,  and  smallest  for  very  rugged  terrain  where  the  contours  are 
very  closely  spaced.  Conversely,  the  problem  of  misregistration  is  most  serious  when  the 
slope  of  the  terrain  is  changing  rapidly  compared  with  the  pixel  size  of  the  Landsat  data. 
In  these  areas  the  method  cannot  be  expected  to  improve  the  accuracy  of  classif ications 
and  indeed  may  even  degrade  them. 

7  CONCLUSIONS 

A  technique  has  been  described  that  can  be  used  to  generate  a  terrain  height  matrix 
from  a  set  of  height  contours.  The  use  of  this  terrain  matrix  to  derive  slope,  aspect, 
and  the  variation  in  direct  solar  illumination  has  also  been  described.  As  an  example 
these  parameters  were  derived  from  data  obtained  from  an  Ordnance  Survey  map  and 
presented  in  the  form  of  grey  tone  images.  Finally,  the  terrain  data  was  used  to 
demonstrate  how  an  attempt  couLd  be  made  to  reduce  the  topographic  effect  on  Landsat 
imagery.  It  was  concluded  that  in  order  to  establish  any  improvement  in  classification 
accuracy  arising  from  the  use  of  terrain  data  with  Landsat  imagery,  a  large  test  area 
is  needed  for  which  accurate  ground  truth  data  is  available.  It  is  expected  that 
attempts  to  use  terrain  data  with  Landsat  and  other  imagery  of  land  surfaces  will  increase 
as  both  the  availability  and  quality  of  suitable  topographic  data  bases  improves  and 
provision  is  made  for  its  routine  use. 


Table  I 


HEIGHT  MATRIX  COMPARED  WITH  SPOT  HEIGHTS  FOR  FIG  1 


Spot  height 
(ft) 

Grid  height 
(ft) 

436.0 

433.7 

343.0 

346.4 

365.0 

359.3 

313.0 

320.6 

353.0 

358.2 

281.0 

294.2 

335.0 

330.7 

318.0 

321.8 

602.0 

600.9 

394.0 

400.4 

379.0 

378.  1 

367.0 

366.4 

442.0 

428.6 

489.0 

481  .7 

414.0 

419.6 

453.0 

436.4 

381  .0 

381  .8 

316.0 

320.9 

377.0 

376.5 

571.0 

566.9 

329.0 

335.7 

Di  £  f erence 
(ft) 


2.3 

-3.4 

5.7 

-7.6 

-5.2 

-13.2 

4.3 
-3.8 

1  .  1 
-6 . 4 


Mean  difference: 


-0.04  ft 


Standard  deviation  of  the  difference:  7.1  ft 
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Table  2 


COMPARISON  OF  CLASSIFICATIONS 


Total  number  of 
pixels  classified 


Data  used 

Original 

Original  ratioed 

Illumination 

corrected 


*  Test  areas  cover  866  pixels. 


Number  of  pixels 
classified  within 
the  test  areas* 


Percentage  of  test 
areas  classified 


Table  3 

PIXEL  STATISTICS  FOR  WOODLAND  AREAS  BEFORE  AND  AFTER  ILLUMINATION  CORRECTION 


D  .  Wavelength  range 

D3nd  /  x 

(um) 


0.5  to  0.6 
0.6  to  0.7 
0.7  to  0.8 
0.8  to  1.1 


Mean 

pixel  intensity 

Variance 

of  pixel  intensities 

Before  correction  After  correction 

80 

79 

17 

137 

55 

54 

24 

60 

90 

88 

109 

105 

102 

99 

245 

151 
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(a)  Height 


(b)  Slope 


(c)  Aspect 


(d)  Variation  in  direct 
solar  illumination 

(The  sun  is  at  an  elevation  of  15° 
and  an  azimuth  of  159°  from  North) 


Fig  2  Grey  tone  image  representation  of  the  four  terrain  parameters  for 
the  central  4  km  square  of  the  test  area  The  pixel  size  is  50  m 
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Landsat  band  7  for  the  central  4  km  square  of  the  test  area  before  and 
correcting  for  the  topographic  effect.  The  pixel  size  is  50  m 
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